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We consider a two-component immiscible Bose-Einstein condensate with dominating intra-species repulsive 
density-density interactions. In the ground-state phase of such a system only one condensates is present. This 
can be viewed as a spontaneous breakdown of Z 2 symmetry. We study the phase diagram of the system at finite 
temperature beyond mean-field approximation. In the absence of rotation, we show that the system undergoes 
a first order phase transition from this ground state to a miscible two-component normal fluid as temperature is 
increased. In the presence of rotation, the system features a competition between vortex-vortex interaction and 
short range density-density interactions. This leads to a rotation-driven “mixing” phase transition in a spatially 
inhomogeneous state with additional broken U(l) symmetry. Thermal fluctuations in this state lead to nematic 
two-component sheets of vortex liquids. At sufficiently strong inter-component interaction, we find that the 
superfluid and Z 2 phase transitions split. This results in the formation of an intermediate state which breaks 
only Z 2 symmetry. It represents two phase separated normal fluids with density imbalance. 


I. INTRODUCTION 

Bose Einstein condensates (BECs) serve as highly useful 
synthetic model systems for a wide variety of real condensed 
matter systems, due to their tunable interactions using mag¬ 
netic and optical Eeshbach-resonances. By creating mixtures 
of the same boson in different hyperfine states, one effec¬ 
tively creates multicomponent condensates [1-4]. Eurther- 
more, by using crossed lasers, one may set up lattice model 
systems with a vast combinations of intersite hopping ma¬ 
trix elements, as well as intrasite interactions, both intra- and 
interspecies [5- 10]. This means that these model systems, 
apart from being interesting in their own right, emulate var¬ 
ious aspects of a plethora of condensed matter systems of 
great current interest, such as multicomponent superconduc¬ 
tors, Mott-insulators, and even topologically nontrivial band 
insulators. The latter follows from the recent realization of 
synthetic spin-orbit couplings in such condensates [11-13]. 
Of particular interest is the physics of these systems in the 
strong coupling regime. 

Spinor condensates with two components of the order pa¬ 
rameter represent a first step away from ordinary single¬ 
component condensates. This extension opens up a whole 
vista of physics which has no counterpart compared to single¬ 
component condensates, due to the wide variety of inter¬ 
species couplings that may be generated. Thus, these syn¬ 
thetic systems display physics which is beyond what is ordi¬ 
narily seen in condensed matter systems. 

The parameter range where inter-component density- 
density interactions exceed intra-component density-density 
interactions signals the onset of immiscibility, or phase sepa¬ 
ration, of the two components. Numerical works solving the 
Gross-Pitaevskii ground-state equations have also found in¬ 
teresting vortex lattices in this regime [14—22]. The effect of 
the repulsive inter-component interactions overpowering the 
intra-component interactions causes the condensate to form 
intertwined sheets of vortices [23]. The condition for immis¬ 
cibility is readily realized experimentally, using magnetic and 


optical Eeschbach resonances [16, 24]. 

In this paper, we present results of large-scale Monte-Carlo 
simulations of a two-component Bose Einstein condensate at 
finite temperature. In a previous work, we have considered 
in detail the effect of thermal fluctuations for the case where 
the inter-component density-density interaction is less than or 
equal to the intra-component interaction [25]. In the present 
paper, we focus on the regime of density-density interactions 
where the inter-component interactions is larger than the intra¬ 
component interactions. This regime is qualitatively different 
from the case in which the intra-component interactions dom¬ 
inate. 

Previous works have studied the effect of an inter¬ 
component density-density interaction on the rotation- 
induced non-homogeneous ground states. These works were 
mostly limited to two spatial dimensions solving the Gross- 
Pitaevskii ground-state equations [14-23]. although certain 
aspects of the three-dimensional case were also studied at 
mean-field level [15, 22]. Here, we extend previous works in 
the immiscible regime to the case of finite temperatures and 
higher dimensions, including the full spectrum of density and 
phase-fluctuation of the condensate ordering fields. 


II. DEFINITIONS 

A. Model 

We consider a general Ginzburg-Landau(GL) model of an 
W-component Bose-Einstein condensate, which in the ther¬ 
modynamic limit is defined as 
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where 
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is the Hamiltonian. Here, the field formally appears as a 
non-fluctuating gauge-field and parametrizes the angular ve¬ 
locity of the system. The fields ipi are dimensionfull complex 
fields, i and j are indices running from 1 to denoting the 
component of the order parameter (a “color” index), ai and 
Qij are Ginzburg-Landau parameters, $o is the coupling con¬ 
stant to the rotation-induced vector potential, and rrii is the 
particle mass of species i. For mixtures consisting of dif¬ 
ferent atoms or different isotopes of one atom, the masses 
will depend on the index i, while for mixtures consisting of 
same atoms in different hyperfine states, the masses are equal 
among the components i. The inter- and intra-component 
coupling parameters gij are related to real inter- and intra¬ 
component scattering lengths, Uij, in the following way 


g^i = 

gij = 
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that (IV'iP) = (|■*/' 2 p) when w > 0. (ai ^ would act as 
an external field conjugate to the pseudo-magnetization of the 
system, and would destroy the Ising-like phase transition we 
report on below). 

Conversely, for a binary mixture of homonuclear cold 
atoms, one may express the ratios of intra- to inter-component 
scattering lengths in terms of w and rj, as follows 

Qi2 ^ l-uj/g 
an 1 + uj/rj 

Note that the ratio of the scattering lengths only depend on the 
ratio oj/rj. 

We discretize the model on a cubic lattice with sides L by 
defining the order parameter field on a discrete set of coordi¬ 
nates —>■ ipr,i, r C (*x + jy + kz\i,j,k = 1,..., L). 

The covariant derivative is replaced by a forward difference. 
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Here, the lattice version of the non-fluctuating gauge field is 
parametrized in Landau gauge, r = (0, 2Trfx, 0), where / 
is the number of vortices per plaquette, or filling fraction. The 
lattice spacing, a, is fixed to be smaller than the characteristic 
length scale of the variations of the order parameter, and ft G 
(x, y, z) is a unit vector. 

Thus, the lattice version of the Hamiltonian we consider is 
given by 


Here, my = mi mj/(mi -f mj) is the reduced mass. In this 
work, we focus on using BECs of homonuclear gases with 
several components in different hyperfine states, i.e. = 
mVi The system we primarily have in mind is a mixture 
of Rb®^ atoms in two different hyperfine states F = 0 and 
F = 1, such that N = 2. Note that when gi 2 = g 2 i = 
^g > gii = 522 = g, i-e. A > l, there is a strong tendency 
in the system to phase separate, leading to two immiscible 
quantum fluids. For a homonuclear binary mixture, such as 
the mixture of Rb®^ atoms mentioned above, we have = 
mi/2. Then, it suffices that > an for the inter-component 
density-density interactions to dominate the intra-component 
density-density interactions. 

In the following, we have introduced dimensionless cou¬ 
pling parameters and fields, following Appendix A of Ref. 25. 
It is convenient to rewrite the potential (repeated indices are 
summed over) 

V = ai + g^j Itpif , (5) 

by introducing interaction parameters rj, uj, such that g = r] + 
w, and Xg = rj—uj, i.e. g = g{l+X)/2, w = ^(l—A)/2. Here, 
A denotes the ratio between the inter- and intra-component 
interactions. Then, Eq. 5 takes the form (up to an additive 
constant) 

V ={ai+ 2p)|V'iP + (a 2 + 2p)|V’2p 

+ viM + - 1 )" + - 1 ^- 2 ^^ ( 6 ) 

For A > 1, a; < 0, with the proviso that ?7 + a; = ?7 — |a;|>0 
for stability. Furthermore, we will assume that ai = a 2 , such 
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Here, we have written the order parameter fields as real am¬ 
plitudes and phases, 1/^,1 = |'0r.i| e*®‘' L In addition, we have 
defined an energy scale, Jq = agO^/^Q, where ag and pg 
are the parameters of the Ginzburg-Landau theory at T = 0. 
Throughout, we fixry = 5.0 and ai + 2g = a2+2,g = 0. This 
guarantees a non-zero ground state condensate density for all 
values of g. 

B. Ground state symmetry 

Eq. 9 defines two superfluids coupled by density-density 
interactions. When there is no phase separation, we have a 
U(l) X U(l) symmetry broken in the ground state. When the 
inter-component interaction is equal to the intra-component 
interaction the system breaks SU{2) symmetry. Here, we are 
interested in the phase separated case. In this case, the sys¬ 
tem breaks an additional Z 2 symmetry, corresponding to in¬ 
terchanging ipi -O' '(/> 2 - That is, when a; > 0, |'0ip = |'i/' 2 p 
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is favored. This represents a Z 2 -symmetric state. On the 
other hand, when uj < 0, |'0ip ^ |'!/' 2 p is favored, such 
that IV'iP — |' 02 p may acquire a nonzero expectation value, 
with equal probabilities that the expectation value is either 
positive or negative. This corresponds to breaking an Ising- 
like Z 2 symmetry. Thus, the ground state breaks a composite 
[/(I) X Z 2 symmetry. 


C. Observables 

The equilibrium phases the model are characterized by sev¬ 
eral order parameters. To identify the Ising-like, phase sepa¬ 
rated order of the system we define 


A = 


(IV'i|^)-(IV'2|^) 


( 10 ) 


To further characterize vortex structures, we examine the 
structure factor of the vortices, defined as 
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This is simply the Fourier-transform of the x-averaged vortic- 
ity. To improve the resolution of the interesting g-vectors, we 
remove the qj^ = 0 point from the hgures. We also compute 
the specific heat capacity 




(17) 


as a means of precisely locating the various transition points. 


where is the thermal and spatial average of 


D. Details of the Monte-Carlo simulations 


= (11) 
r 

A hnite value of A signals relative density depletion in either 
of the condensates. In addition to Z 2 order, it is important to 
monitor the U(l) ordering of the system. The helicity mod¬ 
ulus measures phase coherence along a given direction of the 
system. It is defined as 


1 d^F[0'] 

“ L3 dSf, 

< 5„=0 
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Here, F[9'] is the free energy with an infinitesimal phase twist, 
Sfi, applied along the p-direction, i.e., we make the replace¬ 
ment 


Sr,i 0r,i — 0r,i ~ ^ ' 


(13) 


in F. 

We also identify the nature of the phases by computing 
thermal averages of real-space configurations of densities 
(|'0i(r_L)|^) and vortices (|ni(r_L)|^) in the system. These are 
computed by averaging the quantity along the z-direction of 
the system, with subsequent thermal averaging. That is, 

(n*(r_L)) = (^ X! 

and 

^ z 

The vorticity, ni.r is calculated by traversing a plaquette with 
surface normal in the z-direction, adding the phase differ¬ 
ence 0r+ft,i — 9r,i — on each link. If this plaquette sum 
turns out to have a value outside the primary interval, (—tt, tt], 
2n7r(—2n7r) is added to the sum, which inserts a vortex of 
charge +n{—n) on the plaquette. 


We consider the model on a lattice of size x Ly x L^, 
using the Monte-Carlo algorithm, with a simple restricted 
update scheme of each physical variable, and Metropolis- 
Hastings [26, 27] tests for acceptance. Here, Li is the linear 
extent of the system in the Cartesian direction i G (x, y, z). 
In all our simulation, we have used cubic systems = 
Ly=L^= L, with L G {24,32,40,48, 56,64,96,128}. At 
each inverse temperature, 10® Monte-Carlo steps are typically 
used, while 10® additional sweeps are used for equilibration. 
Each Monte-Carlo step consists of an attempt to update each 
amplitude and phase separately in succession, at each lattice 
site. To improve acceptance rates, we only allow each update 
to change a variable within a limited interval around the previ¬ 
ous value, the size of which is chosen by approximately maxi¬ 
mizing acceptance rates and minimizing autocorrelations. The 
Mersenne-Twister algorithm is used to generate the pseudo¬ 
random numbers needed[28]. To ensure that the state is prop¬ 
erly equilibrated, time series of the internal energy measured 
during equilibration are examined for convergence. To avoid 
metastable states, we make sure that several simulations with 
identical parameters and different initial seeds of the random 
number generator anneals to the same state. Measurements 
are post-processed with multiple histogram reweighting[29]. 
Error estimates are determined by the jackknife method[30]. 


III. PHASE DIAGRAM IN THE ABSENCE OE ROTATION 

The model has aU(l) x U(l) x Z 2 symmetry and hence 
the possibility for several phase transitions. However, for the 
parameter set which was simulated a single first-order phase 
transition is found, as we illustrate below for w = —2. Eor 
/ = 0, we have considered L G (24, 32,40,48, 56, 64}, and 
we present the results for L — 64. 

Eig. 1 shows the specific heat, Ising order parameter A, 
and helicity moduli (phase-stiffness or equivalently superfluid 
density) of both components for L = 64 and w = — 2. The 
onset of the Ising-like order parameter A is located at the 
same temperature as the 5-function anomaly of the specific 
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FIG. 1; The specific heat. A, and helicity moduli ^ for 
L = 64 for w = —2. The specific heat features a i5-function 
anomaly, characteristic of a first-order phase transition. A 
has an onset at the same value that the specific heat anomaly 
is found. The helicity modulus of the condensate with the 
smallest density is zero for all values of /3. Note also the 
discontinuities in the Z 2 -order parameter A, and the U(l) 
-order parameters i. The origin of the first-order character 
of the composite Z 2 x U(l) phase-transition is due to 
interaction between domain walls and vortices. Details are 
explained in the text. 


heat. Note that A vanishes discontinuously as f3 approaches 
the transition point from above. Fig. 1 also shows the helicity 
moduli of both components above and below the Ising tran¬ 
sition. In the high-temperature Z 2 -symmetric phase, where 
both components have equal densities ^ = 0. 

The onset of ^ signals a broken U(l) symmetry. Systems 
simultaneously breaking U{1) and Z 2 symmetries, can have 
two independent phase transitions. These are driven by pro¬ 
liferation of vortex loops for the U(l) transition and prolifer¬ 
ation of domain walls for the Z 2 transitions. For the system in 
question, these transitions are not independent. Domain wall 
excitations interact with vortices, and therefore proliferation 
of these topological defects are not independent processes. In 
our system the non-trivial interplay between the U(l) - and 
Z 2 - sectors leads to a single phase-transition. A different ex¬ 
ample of system where interacting Z 2 and U(l) sectors lead 
to a nontrivial phase diagram, is the case of phase-frustrated 
multiband superconductors [31, 32]. 

Consider lowering the temperature from the fully symmet¬ 
ric phase where A = ^ = 0. The Z 2 -symmetry is bro¬ 

ken at a certain temperature such that A ^ 0, i.e. (|V’i|^) 
(|'(/> 2 |^)- Thus, one component gets a reduced average density 
and one gets an enhanced average density. These densities 
determine the bare phase-stiffnesses in the problem. Thus, the 
component with the largest density effectively has less phase 
fluctuations than the other one. Furthermore, due to suppres¬ 
sion of one of the components, the helicity modulus belonging 
to the dominant component becomes non-zero, while the he¬ 


licity modulus belonging to the minor component can remain 
zero. This effect is enhanced as temperature is lowered, since 
A increases monotonically with /3, thus decreasing bare phase 
stiffness of the minor component as /3 is increased. The low- 
temperature phase is therefore a two-component Bose fluid 
where one component is in a U(l) -ordered phase-coherent 
state and the other is in a U(l) -disordered phase-incoherent 
state. The spontaneous appearance of a disparity in densities 
reflects a broken Z 2 symmetry. 

Breaking a Z 2 - or U(l) -symmetry is usually associated 
with a second order phase transition, i.e. a critical point. 
We next discuss how the above situation instead leads to a 
first-order phase transition. Consider heating up the system 
from deep within the low-temperature phase, described in the 
previous paragraph, until the system approaches the vicinity 
of the transition to the fully symmetric phase. Formation of 
domain walls implies regions of suppressed density as well 
as imposes cutoffs to vorticity, thus enhancing phase fluctu¬ 
ations. On the other hand, thermally-induced vortices can 
decrease domain-wall tension. Thus, both the Z 2 -transition 
and U(l) transitions can take place preemptively, compared 
to single-component models. Under such circumstances two 
phase transitions can merge into a single first order phase tran¬ 
sition. 

A similar mechanism for producing a single first-order 
phase transition by merging second-order phase transitions, 
has previously been discussed in different systems with com¬ 
peting topological excitations, such as the so-called s + is su¬ 
perconducting states [31], in U{1) x U{1) systems with elec¬ 
tromagnetic and drag couplings[33, 34], and as an interpre¬ 
tation of observed first order phase transitions in multicom¬ 
ponent gauge theories [35-38]. The phenomenon therefore 
appears quite generic in the physics of multi-component con¬ 
densates. 

To corroborate the above discussion, we have performed 
finite size scaling of the specific heat peak height for system 
sizes L G {24,32,40,48,56,64}. The scaling is shown in 
Fig. 2 on a logarithmic scale, and the exponent obtained is 
ajv = 3.38(8). This is consistent with a first order transition, 
where the exponent would be 3. 

For L — 64, we have performed similar computations for 
w = (-0.1,-0.25,-0.5,-1,-2}. The results for / = 
0 are summarized in Fig. 3. The low-temperature phase 
is a two-component immiscible (phase-separated) superfluid 
with spontaneously broken Z 2 -symmetry, while the high- 
temperature phase is a two-component miscible normal fluid. 
The solid line separating these two phases is a first-order 
phase transition where a U(l) - and a Z 2 -symmetry are bro¬ 
ken simultaneously. The endpoint of the phase-transition line 
at a; = 0 is a point where the system acquires an SU (2)- 
symmetry (see e.g.[25]). 


IV. MIXING AND SUPERFLUID PHASE TRANSITIONS IN 
THE PRESENCE OF ROTATION 

We next consider the effect of imposing a finite rotation on 
the condensate. Our main results are presented for a system 
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FIG. 2: Log-Log plot of the finite size scaling of the height of 
the specific heat curves at L S {24,32,48,56,64}. The 
parameters are w = — 2.0, / = 0, and rj = 5.0. The exponent 
obtained is a = 3.38(8). 
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FIG. 3: The phase diagram of the two-component Bose 
Einstein condensate at zero rotation, / = 0, p = 5, and 
w < 0. The high-temperature phase, region I, is a 
two-component miscible normal fluid. The low-temperature 
phase, region II, is an immiscible (phase-separated) 
superfluid. The solid line separating these two phases is a 
first-order phase transition where a composite U(l) x Z 2 
-symmetry is broken. 


size of L = 64 and / = 1/32, but we have considered system 
sized L e {32, 64, 96,128}. 

Introducing a finite amount of (rotation-induced) vortices 
in the ground state significantly alters the simple phase di¬ 
agram presented in Fig. 3. The first effect is to suppress 
the parameter regime where a broken Z 2 -symmetry is found, 
A ^ 0. Recall that for f — 0, any w < 0 sufficed to bring 
about A 7 ^ 0 at sufficiently low /?, as seen in Fig. 3. A fi¬ 
nite amount of vortices alters this. Vortices interact via long 


range current-current interactions. It is energetically favor¬ 
able to maximize the distance between vortices, subject to 
the constraint that a specific number of them has to be con¬ 
tained within a given area perpendicular to the direction of 
rotation. This effect leads to a uniform distribution of min¬ 
ima (equivalently maxima) in the condensate densities. On 
the other hand, density suppression by vortices in one com¬ 
ponent in general allows the second to nucleate. The short- 
range repulsive inter-component density-density interaction 
(ry — a;)(|'!/'ip|'!/ 2 p + |'*/' 2 p|' 0 iP) (which exceeds the intra¬ 
component density-density interaction (77 + -f 

IV' 2 p|' 02 p) for uj < 0), tends to produce regions where den¬ 
sities one component is large while the other is small, and 
vice versa. Below a critical value of —w = —tUc « -1-0.6, 
we do not see any onset of A 7 ^ 0 at any value of /3 as the 
system is cooled from a uniform state. That is, the interface 
tension between the phases is sufficiently low and the overall 
free energy, which includes long range inter vortex interaction 
is minimized by the state with A = 0. 

For the subsequent discussion, it helps to consider a 
schematic phase diagram of the system with / 7 ^ 0 , which we 
have obtained through large-scale Monte-Carlo simulations. 
The phase diagram is shown in Fig. 4. Region I denotes 
the simple translationally invariant high-temperature Z 2 - and 
U(l) X U(l) -symmetric two-component phase with equal 
densities of both condensate components. Region II shows 
the Z 2 -symmetric striped phase. Region III is a region with 
broken Z 2 -symmetry, with one stiff condensate component in 
a uniform hexagonal vortex lattice phase, and one component 
in a uniform vortex liquid phase. Region IV is a region with 
broken Z 2 -symmetry, but with the two condensates both in 
a vortex-liquid phase. Thus, the phase transition separating 
region I from region II is a phase-transition line separating a 
two-component isotropic vortex liquid from a two-component 
striped (nematic) vortex liquid. The line separating region I 
from Region IV is one where a Z 2 -symmetry is broken, and 
the line separating region II from region IV is one where a 
translational symmetry is broken and the system acquires non¬ 
zero helicity modulus 


A. Transition from region I to region II 

We first consider the thermally driven transition from 
the high-temperature symmetric two-component vortex liq¬ 
uid phase, region I, to the low-temperature two-component 
striped (nematic) phase, region II, for fixed negative to, but 
where |w| < \ujc\, i.e. to the left of the splitting point where 
Region IV opens up. 

In Fig. 5 we show the specific heat cy, helicity moduli in 
the z-direction ^ as the inverse temperature /3 is varied, for 
/ = 1/32 and uj — —0.50. This corresponds to a value of 
—u! to the left of the splitting point where Region IV opens 
up (see Fig. 4). The longitudinal helicity moduli ^ of both 
components develop a finite expectation value. The onset of 
this finite value is accompanied by an anomaly in the specific 
heat. 

We note the sharp, (5-function anomaly in the specific heat 
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FIG. 4; The phase diagram of the two-component Bose 
Einstein condensate at finite rotation, f ^ 0, rj = 5, and 
w < 0. Negative oj may lead to the breaking of the Z 2 
-symmetry in the problem, in addition to the usual breaking 
of the obvious U(l) x U(l) - symmetry. Region I is a Z 2 - as 
well asU(l) X U(l) -symmetric two-component 
vortex-liquid phase. Region II is a Z 2 -symmetric striped 
(nematic) phase consisting of a two-component vortex liquid 
with broken translational symmetry in a direction 
perpendicular to the stripes, but not in the direction parallel 
to the stripes. Region III is a phase with broken Z 2 
-symmetry, and with broken translational symmetry in one 
condensate component, but not the other. Region IV is 
similar to Region III, except that no translational symmetry is 
broken in either condensate component. Details are 
explained in the main body of the paper. 
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FIG. 5: Specific heat, cy, and helicity moduli along the 
z-axis ,Tz i, with / = 1/32 and w = —0.50, i.e as the 
system transitions from Region I to Region II in Fig. 4. 


FIG. 6; Histograms of the probability distribution of the 
internal energy per site, U/L^, at the transition point 
uj = —0.5, P « 0.9995, separating Region I from Region II 
in Fig. 4, for L = {64,96}. multi-histogram reweighting was 
used to obtain histograms with approximately equal peak 
heights. 


and the discontinuous behavior of the helicity moduli in both 
components. These features are all straightforwardly inter¬ 
preted as signals of a first-order phase transition. This is fur¬ 
thermore borne out by performing a computation of the his¬ 
togram of the free energy versus internal energy of the system 
at precisely at the transition, see Fig. 6. This shows a double¬ 
dip structure with a peak in between, the standard hallmark of 
two degenerate coexisting states separated by a surface whose 
energy is given by the height of the peak between the min¬ 
ima. This surface energy clearly scales up with system size 
(more precisely it scales with the cross-sectional area of the 
system), while the difference between the energies separating 
the two degenerate states approaches as finite value as sys¬ 
tem size increases. The histograms develop into two separate 
(5-function peaks as the system size increases, while the differ¬ 
ence in the internal energy between the two degenerate states 
of equal probability (equivalently of equal free energy) is the 
latent heat of the system. The latter clearly approaches a h- 
nite value per degree of freedom as the system size increases, 
demonstrating the hrst-order character of the transition. 

To further characterize the transition I —> II, Fig. 7, 
shows the Z 2 order parameter A and the structure func¬ 
tions G (1,2) in a narrow range around the tran¬ 

sition point. From the top panel, it is seen that A = 0 
for all P considered. Moreover, we see that as /3 is in¬ 
creased, the structure function evaluated at a q-vectors, Gc = 
(±7r/32, ±7r/32) on the Bragg-circle of the vortex-liquids 
are reduced, while the structure function evaluated at Bragg 
peaks, Gs = (±7r/32, =F7r/32), corresponding to a uniform 
striped phase of / = 1/32 increases. The onset of the lat¬ 
ter marks the transition from a uniform two-component vor¬ 
tex liquid to a two-component nematic vortex liquid, a striped 
phase. The mechanism for producing the striped phase is de- 
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scribed above. Note that in the thermodynamic limit, isolated 
vortex sheets can be expected to be in the state of one dimen¬ 
sional liquid at any finite temperature in analogy with the ab¬ 
sence of crystalline order in one dimensional systems. 

We thus conclude that the transition from Region I to re¬ 
gion II is a first order phase-transition involving the break¬ 
ing of a composite U(l) x U(l) symmetry, from an isotropic 
two-component vortex liquid in Region I to a two-component 
nematic phase of intercalated lattices of stripes of one¬ 
dimensional vortex liquids in Region II. We next go on to con¬ 
sider in some more detail the structure functions, primarily to 
gain more insight into the character of the striped phase of 
Region II. 

The four bottom panels of Fig. 7 show the structure func¬ 
tions e (1,2) at two values of /3, /3 = 0.990 and 

P = 1.010. At /3 = 0.990, both structure functions show 
ring-like structures characteristic of an isotropic liquid. No¬ 
tice also that the intensity of the rings are equal, which is a 
consequence of the fact that A = 0. At /3 = 1.010, both 
structure functions have developed Bragg peaks in one direc¬ 
tion bot no Bragg peaks in the corresponding perpendicular 
direction. This is indicative of a striped phase. 

This may be further corroborated by correlating the struc¬ 
ture functions with real-space vortex structures for various 
values of /3. This is shown in Fig. 8 

One aspect of the structure functions shown in the two bot¬ 
tom rows of Fig. 8, is particularly important. Consider first 
the case /? = 0.900, well within region 1 for w < Wc. This 
is shown in the hrst column of Fig. 8. The real-space vortex 
configurations in both components are disordered. Moreover, 
'S'i(q_L),* € (1,2) both feature ring-structures characteristic 
of an isotropic liquid phase. The value of |q| at which the 
rings appear is a measure of the average inverse separation 
between the vortices in the isotropic liquids. The intensities 
of both structure functions is the same. Consider next the case 
P = 0.995, shown in the second column of Fig. 8. From the 
real-space pictures, one discerns a tendency towards stripe- 
formation. This is reflected in ^^(qj,),^ S (1,2), where the 
ring-like structures now instead are anisotropic, developing 
peaks in the direction perpendicular to the direction of the in¬ 
cipient stripes. At even lower temperatures p = 1.100, well 
within region II where stripes are fully developed, the ten¬ 
dency towards anisotropies in S'i(q_L), i G (1, 2) is even more 
obvious. This is shown in the last column of Fig. 8. In this 
case, Bragg-peaks have fully developed in the directions per¬ 
pendicular to the stripes. There are, however, no Bragg peaks 
in the direction parallel to the stripes. If the stripes were per¬ 
fectly straight, there would be two weak Bragg peaks in these 
directions. This would be the one-dimensional analog of the 
ring-like liquid structures of the isotropic liquid. The value of 
|q| at which this single weak peak occurs corresponds to the 
inverse average separation between vortices within the stripes. 
The reason they are not observed in our calculations, is due to 
the slight fluctuations in the shape of the stripes, which wash 
the Bragg peaks out. 

We thus conclude that region II is a striped phase where 
the stripes form one-dimensional (ID) vortex liquids. Vor¬ 
tices in quasi-ID systems have finite energy and cannot form 
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FIG. 7: Z 2 order parameter A and vortex structure functions 
‘S'i(qj.), * G (1, 2) in the vicinity of the transition from 
Region I to Region II, Fig. 4, with / = 1/32, w = —0.50. 
Panel (a) shows A as a function of P, as well as structure 
functions at specific points in reciprocal space, Si{Gs) and 
5'i(Gc). The four bottom panels show the structure functions 
S'i(qj^) and S' 2 (q_L) for the two values P = 0.990 and 
P = 1.010. 


a ID solid at any finite temperature. This is consistent with 
the structure factor we observe. On the other hand, the inter¬ 
action between stripes may not be negligible, so the details of 
the phase diagram in Region II warrant further investigation. 
A notable feature of this state is the finite helicity moduli in z- 
direction, even if the structure factors show absence of vortex 
ordering within stripes. This highly unusual situation origi¬ 
nates with the positive interface energy between the two con¬ 
densates. That is, consider a stripe-liquid in x-direction. A 
vortex line in the z-direction is free to execute transverse me- 
anderings in the x-direction. A superflow in the z-direction 
would produce a y-component of the Magnus-force on the x- 
components of the fluctuating vortex lines. However, vortex 
segments are restrained from moving in y-direction due to the 
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inverse temperature is varied between each column, (3 G {0.900, 0.995,1.100}. Each row show, in order, averaged vortex 
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stripe interface tension. This results in the observed finite he- 
licity modulus in z-direction. Similar results are found for a 
number of other w-values we have considered, for —oj < 0 . 6 . 
[39] 


B. Transition from region I to region III, via region IV 

Increasing —w further, such that the inter-species density- 
density interaction increases, eventually favors a different pat¬ 
tern of phase-separation of the two components, despite the 
uniforming effect of long-range current-current interactions 
between rotation-induced vortices. This leads to a broken 
Z 2 -symmetry. The condensate component with a globally 
suppressed density will therefore be in a vortex-liquid phase 
while the condensate component with globally enhanced den¬ 
sity will be in a vortex lattice phase. The combined preemp¬ 
tive U(l) X Z 2 phase transition found for / = 0, now splits 
into two separate phase transitions. The splitting occurs be¬ 
cause the U(l) -sector directly couples to the rotation, while 
the Z 2 -sector does not. The phase-transition in the stiff U(l) - 
sector, which is a vortex-lattice melting, is therefore separated 
from the Z 2 -transition by an amount which depends on /. 

Since A increases with —w beyond —uJc, the transition 
temperature for the Z 2 -transition increases. This effectively 
makes the dominant component stiffer as —uj increases. Thus, 
the temperature for melting the vortex lattice in the dominant 
U(l) -sector also increases with increasing —w. The splitting 
between the U(l) - and the Z 2 -sectors also increases as —w 
increases. 

This is illustrated in Fig. 9, showing A, specific heat Cy, 
and 1 , T ^_2 as functions of /3. The Z 2 order parameter A 
has an onset at , at which the specific heat has an anomaly. 
There is no onset of 1 , showing that component 1 remains 
in a vortex liquid phase. Component 2 forms a vortex solid 
at lower temperature, as evidenced by the onset of 2 - This 
happens at a /3u(i) which is separated from as explained 
above. 

Fig. 10 shows the structure functions 5'i(q_L) and S' 2 (qj_) 
at u! = —4.0 at three different values of (3, namely /3 = 
(0.275, 0.290, 0.310). These values correspond to Regions I, 
IV, and III in Fig. 4, respectively. Here again, we see the 
freezing of one component across the transition, while the 
other component remains in the liquid phase. The additional 
information we get out of these panels is that one compo¬ 
nent remains an isotropic vortex liquid, while the other com¬ 
ponent freezes into a hexagonal vortex liquid. This sets the 
low-temperature Region III (see Fig. 4) at w = —4.0 drasti¬ 
cally apart from the low-temperature Region II (see Fig. 4) 
at OJ = —0.50. The latter features a low-temperature two- 
component nematic vortex liquid phase with broken rotational 
invariance, the former case features a low-temperature mixed 
isotropic vortex liquid/hexagonal vortex lattice phase. 

For a more detailed overview of the transition. Fig. 11 show 
the evolution of {ni{r±)), (|' 0 i(rj_)|), and S'i(qj^) across the 
three regions, I, IV, and III. If one follows the evolution of the 
vortex densities in each component, it is seen that the com¬ 
ponent which acquires a low stiffness in region IV and HI is 



/3(l/Jo) 


FIG. 9; The phase transitions between Region I and Region 
IV, and between Region IV and Region III, for 
/ = 1/32, w = —4.0, and L = 128. Note the separation 
between the onset of A and T^ i. The onset of A signals the 
breaking of a Z 2 -symmetry, along with the associated 
anomaly in specific heat Cy- This marks the transition from 
Region I to Region IV in Fig. 4. In Region IV, we have 
A ^ 0, while both components remain in isotropic vortex 
liquid states. In passing from Region IV to Region III in Fig. 
4, the onset of one of the helicity moduli, say, signals 
the freezing of the vortex liquid in the corresponding 
component, while the absence of an onset of the helicity 
modulus, 2 say, in the other component shows that this 
component remains in a vortex liquid phase. The onset of 
T 2 1 signals the breaking of a U(l) -symmetry associated 
with vortex-liquid freezing. 


virtually unchanged, i.e. it remains in a completely uniform 
state. The other component, on the other hand, evolves from a 
uniform state in region I, through being close to freezing into 
a hexagonal lattice in region IV, and finally into a hexagonal 
structure in region III. The amplitude densities corroborate 
this picture. In region I they are on average equal and uni¬ 
form, while in region IV the difference in stiffness is clearly 
seen. Here some inhomogeneities arise in the stiff component 
as the vortices are close to entering a hexagonal phase, which 
is also reflected in the soft component simply because of the 
local intercomponent repulsion. In region III, the amplitude 
density of the stiff component is high and uniform with small 
dips corresponding to each vortex. The soft component is low 
and uniform with small peaks, again due to intercomponent 
interactions. 


C. Transition from region II to region III 

Finally, we consider the transition from Region II to Re¬ 
gion III. In Region II, we have A = 0, while in Region III, 
A 7 ^ 0. Therefore the Regions II and III are separated by a Z 2 
-symmetry breaking. Stripe-forming systems in general have 
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FIG. 10; The phase transitions of the system for 
/ = 1/32,0; = —4.0. Structure functions 5'i(qj_) and 
S' 2 (q_L) at three different values of P, namely 
(3 = (0.275,0.290, 0.310), corresponding to Regions I, IV, 
and III in Fig. 4, respectively. 


complicated structural transitions. We find an intermediate 
regime where the lattice of stripes has disordered, but where 
the hexagonal lattice/isotropic liquid-mixture has not yet fully 
developed. This results in multiple metastable, but robust co¬ 
existing phases of vortices in components 1 and 2 residing in 
different parts of the condensate. These two coexisting phases 
are separated by a surface of positive surface energy. This sur¬ 
face constrains the motion of vortex systems. As a result, in 
the finite systems which we simulate, the helicity moduli 
acquire nonzero values in both components in this intermedi¬ 
ate regime. 

As —uj is increased further, such that one component be¬ 
comes dominant and the other is suppressed, the minor com¬ 
ponent becomes normal. Note that when the inclusions of 
the normal component become isolated, they represent quasi- 
ID subsystems. Quasi-ID systems are superfluid only at zero 
temperature. However, simulations on finite systems may still 
display finite helicity modulus. As the density of the com¬ 
ponent increases, the corresponding intra-component current- 
current interaction between the rotation-induced vortices in 
this component increases. Hence, the intra-component long- 
range interaction for this component dominates, and a hexag¬ 
onal vortex-lattice results. Consequently, the helicity-moduli 
in the two components have quite different behavior as —uj 
increases. In the component that eventually takes up a vortex 


lattice state, it increases monotonically with —u. In the other 
component, it is non-monotonic as a function of w, eventually 
approaching 0 deep into Region III. 

Typical examples of the vortex structures that appear be¬ 
tween Region II and Region III in Fig. 4 are shown in Fig. 12. 
These are all metastable, long-lived states which prevent equi¬ 
libration of the system. We have been unable to locate the 
sharp separatrix between these two regions, and whether there 
are other stable intermediate phases due to the lack of equi¬ 
libration. Note that this problem is known in other stripe¬ 
forming systems where phases are separated by metastable 
and glassy states [40, 41], 


V. CONCLUSIONS 

In this paper, we have considered the states of a two- 
component Bose Einstein condensate in the situation where 
inter-component density-density interactions dominate over 
the intra-component density-density interactions. The two 
components of the condensate are assumed to be comprised 
of homonuclear atoms in two different hyperfine states. The 
problem features an Ising-like symmetry. This Ising (or 
Z 2 ) symmetry emerges from the dominance of the inter¬ 
component interactions over the intra-component ones. The 
spontaneous breaking of this Ising-symmetry corresponds to 
a spontaneously generated, interaction-driven, imbalance be¬ 
tween condensates in different hyperfine states. 

At finite rotation, we find four regions, denoted Regions I, 
II, III, and IV, of thermodynamically stable states, see Fig. 4. 
Region I is a high-temperature regime where the system re¬ 
mains in a two-component isotropic vortex liquid phase with 
equal densities of both components, i.e. no imbalance be¬ 
tween condensate components in different hyperfine states. 
Region II is a nematic phase (broken rotational symmetry) 
with ordered stripes of one-dimensional vortex liquids, and 
with no imbalance between components of different hyper¬ 
fine states. This state features a spontaneously broken com¬ 
posite U(l) X U(l) -symmetry, but is Z 2 -symmetric. In 
addition it spontaneously breaks translation symmetry in one 
direction due to formation of periodic modulation of conden¬ 
sates. Region III is a mixed state with one component in a 
U(l) -symmetric isotropic vortex liquid phase while the other 
component resides in a hexagonal vortex lattice phase with 
broken U(l) -symmetry. The origin of the different behaviors 
of the two components is that Region III also features a spon¬ 
taneously broken Z 2 -symmetry, i.e. an imbalance between 
the density of one hyperfine state and the other. The compo¬ 
nent with a large density has higher phase stiffness than the 
component with the lower density, hence the discrepancy in 
their vortex states. Finally, Region IV is a region intermediate 
between Region I and Region III, in which U(l) -symmetry 
is not broken in either of the components, but where a spon¬ 
taneously generated imbalance between densities of hyperfine 
states exists. Both components are in an isotropic and disor¬ 
dered vortex state. 

The phase transition from Region I to Region II in Fig. 
4 is a first-order composite U(l) x U(l) transition. The 
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FIG. 11: Tableaux showing detailed real space pictures of the transition from region I to region III, via Region IV. The inverse 
temperature is varied between each column, (3 G {0.275,0.290, 0.310}, and w = —4. Each row shows, from top to bottom, 
averaged vortex densities of each component, (ni(r±)), and averaged amplitude densities of each component, (|t/;i(rj^)|). Note 
the vortex-ordering in one of the components, and the lack of vortex-ordering in the other component, as the system transitions 
from the symmetric phase region I (/3 = 0.275) to the low-temperature phase region III = 0.310). Note also the disparity in 
density-amplitudes in the two components in the intermediate regime region IV (/3 = 0.290), due to the Z 2 -symmetry breaking. 


phase transition between Region I and Region IV is asso¬ 
ciated with a spontaneous Z 2 symmetry breaking where a 
density-imbalance between condensates of different hyper- 
fine states sets in. The phase-transition between Region IV 
and Region III is a first order U(l) transition associated 
with the freezing of an isotropic vortex liquid in one compo¬ 
nent into a hexagonal vortex lattice in the same component, 
while the other component (the one with depleted density 
due to the Z 2 -symmetry breaking) remains in the isotropic 
vortex liquid phase. The phase transition from Region II 
to Region III, driven by increasing the dominance of inter¬ 
component density-density interactions over intra-component 
density-density interactions, involves at the very least a spon¬ 
taneous breaking of a Z 2 -symmetry as the two condensate 
components pass from a nematic state of intercalated lattices 


of one dimensional vortex liquids into a mixed state of an 
isotropic vortex liquid in one component and a hexagonal vor¬ 
tex lattice in the other component. Other than that, this tran¬ 
sition is characterized by a broad regime of metastable states 
with inhomogeneous phase separation. 

Fig. 11 suggests that the rotation frequency is much smaller 
than the second critical frequency flc 2 - A very rough esiti- 
mate, based on core size, gives cx 0.1Oc2. This puts the sys¬ 
tem well outside the regime of lowest-Landau level physics. 
The system is therefore indeed in a regime where it makes 
sense to talk about vortex-degrees of freedom rather than ze¬ 
roes of the order parameter as the relevant degrees of freedom. 
For this rotation frequency, we have found the critical value of 
uj (one of our interaction parameters) to observe phase IV to 
be Wc « —0.6. From this, we find from Eq. 7 that this re- 
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FIG. 12: Tableaux showing detailed real space and reciprocal space pictures of the transition from region II to region III. The 
parameter uj is varied between each column, w S {—0.300, —0.500, —0.700}. Each row shows, from top to bottom, averaged 
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quires scattering lengths 012/011 > 1-3. Since these scatter¬ 
ing lengths a priori are very similar, and can be manipulated 
with Feshbach resonances, it seems feasible to be able to ob¬ 
serve phase IV. In order to see the striped ground states phase 
II, the requirement is only that 012/011 > 1, which certainly 
seems to be within the realms of possibility. 
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